# Loading data

# experimental data
primary <- read.csv("../raw_data/primary_screen.csv")
speciesonly <- read.csv("../raw_data/species_only_data.csv")
indiv <- read.csv("../raw_data/individual_data.csv")

# re-create column for individualID
indiv <- indiv %>%                              # Create numbering variable
      group_by(PaperID) %>%
        mutate(IndividualID = paste(row_number(),unique(PaperID), sep="_"))

# host data
tree <- read.nexus("../raw_data/upham_tree_666.nex")
tax <- read.csv("../raw_data/taxonomy_mamPhy_5911species.csv")

bats <- tax[tax$ord=="CHIROPTERA",]
rodents <- tax[tax$ord=="RODENTIA",]

bat_tree <- drop.tip(tree, setdiff(tree$tip.label, bats$Species_Name))
rodent_tree <- drop.tip(tree, setdiff(tree$tip.label, rodents$Species_Name))

hosts <- read.csv("../raw_data/HostNames_NCBI_Upham.csv")
hosts$Host_Upham <- gsub(" ", "_", hosts$Host_Upham)

# virus names and taxonomy
vtax <- read.csv("../clean_data/Virus_taxonomy.csv")
viruses <- read.csv("../raw_data/VirusNames_translation_Feb23_2024.csv")
viruses <- left_join(viruses, vtax)

#subset to only host reported and upham
hosts <- hosts %>% select(Host_reported, Host_Upham) %>% unique()

# add names
primary <- dplyr::left_join(primary, viruses)
primary <- dplyr::left_join(primary, hosts)
speciesonly <- dplyr::left_join(speciesonly, viruses)
speciesonly <- dplyr::left_join(speciesonly, hosts)
indiv <- dplyr::left_join(indiv, viruses)
indiv <- dplyr::left_join(indiv, hosts)

# merge speciesonly with individual data
dat <- full_join(indiv, speciesonly)

# removing hosts and viruses with NA names
dat <- dat[!is.na(dat$Virus_ICTV),]
dat <- dat[!is.na(dat$Host_Upham),]

# make subset tree to sampled hosts
host_tree <- drop.tip(tree, setdiff(tree$tip.label, dat$Host_Upham))


# These are models of disease, so remove studies where pathology was not reported
dat <- dplyr::left_join(dat, unique(select(primary, PaperID, Virus_reported, Host_reported, PathologyReported_YN)))
# dim(dat)# 2923
dat <- dat[dat$PathologyReported_YN=="Y",]
# dim(dat)# 2726

# subset to only suscptible individuals / species
dat <- dat[dat$Susceptible_YN=="Y",]
# dim(dat) # 2238

# make virus tree
vtax <- as.data.frame(unclass(vtax), stringsAsFactors=TRUE)
frm <- ~superkingdom/realm/kingdom/phylum/class/order/family/genus/Virus_ICTV
vtree <- as.phylo(frm, data = vtax, collapse=FALSE)
vtree$edge.length <- rep(1, nrow(vtree$edge))

# include only viruses in subset data (e.g. after removing mole)
vtree <- drop.tip(vtree, setdiff(vtree$tip.label, dat$Virus_ICTV))
# plot(vtree)


# add host weight and cleaned dose and inoculation route data
dose_dat <- read.csv("../clean_data/dose_data.csv")
dat <- select(dat, -c(Dose_amount, Dose_unit))
dat <- dplyr::left_join(dat, dose_dat)


# adding host specificity data
host_range <- read.csv("../clean_data/host_range_data.csv")

dat <- dplyr::left_join(dat, host_range)
# dim(dat)# 2238

Priors

Visualization of brms default priors (blue) compared to custom priors (yellow) inspired by McElreath’s book Statistical Rethinking.

Priors for the slope parameters:

inv_logit <- function(x) {exp(x)/(1+exp(x))}

prior <- runif(n=10000, min=-100, max=100) # brms default is actually -Inf to +Inf
p <- inv_logit(prior)
prior2 <- rnorm(n=10000, 0, 1.5) # McElreath
p2 <- inv_logit(prior2)

plot(density(p, adj=0.1), xlim=c(0,1), col="#56B4E9", main="", xlab=expression(paste("inv_logit (",beta,")")))
lines(density(p2, adj=0.1), col="#E69F00")

The default improper uniform prior used by brms is flat on the log-odds scale, so when converted to the probabilty scale, puts excess mass near 0 and 1. A Normal(0, 1.5) prior is more reasonable, having a flatter distribution with a mode at 0.5 on the probablity scale.

Priors for the family-level variance parameters:

prior <- rstudent_t(n=10000, df=3, mu=0, sigma=10) # brms default
prior <- prior[prior>0]
# p <- inv_logit(prior)

prior2 <- rnorm(n=10000, 0, 1) # McElreath
prior2 <- prior2[prior2>0]
# p2 <- inv_logit(prior2)

plot(density(prior2, adj=0.1), xlim=c(0,10), col="#E69F00", main="", xlab=expression(paste(sigma)))
lines(density(prior, adj=0.1), col="#56B4E9")

The brms prior is extremely long tailed. Considering we are fitting a non-linear model in which changes in log-odds over over 4 have a diminishing influence due to the ceiling effect of the binomial distribution. With this in mind, a prior of Normal(0,1) constrains the posterior to a more realistic range, and is likely to improve sampling efficiency.

Species level models

# create a host-virus-study level dataset
spec <- dat
spec <- spec[spec$Susceptible_YN==1,]
spec <- spec[!is.na(spec$Susceptible_YN),]


spec %>% group_by(PaperID, Virus_ICTV, Host_Upham, Host_order, Reservoir_match, ses_pd, EvoIso) %>% summarise(n_indivID=n_distinct(IndividualID), N_individuals=sum(unique(N_individuals), na.rm=TRUE)) %>% View()

# any disease, maximum severity, any mortality, number of susceptible individuals
spec <- spec %>% group_by(PaperID, Virus_ICTV, Host_Upham, Host_order, Reservoir_match, ses_pd, EvoIso) %>%
                  summarise(n_individuals = sum(
                            c(n_distinct(IndividualID, na.rm=TRUE), unique(N_individuals)), na.rm=TRUE),
                            n_studies = n_distinct(PaperID), 
                            Disease_YN=max(Disease_YN, na.rm=TRUE),
                            severity_rank=max(severity_rank, na.rm=TRUE),
                            Mortality_YN=max(Mortality_YN, na.rm=TRUE)) %>% unique()

# for combinations in which mortality was all NA max() returns -Inf
# change this back to NA
spec$Mortality_YN[is.infinite(spec$Mortality_YN)] <- NA

# table(spec$Disease_YN)
# 123/(70+123) # 63.7% of h-v-s combinations reported disease

range(spec$n_individuals, na.rm=TRUE)
# 1 - 492


# additional columns for modelling separate phylo and non-phylo species-level effects
spec$Host_name <- spec$Host_Upham
spec$Virus_name <- spec$Virus_ICTV


# scaling / transform continuous  predictors
spec$ses_pd <- scale_half(spec$ses_pd)
spec$EvoIso <- scale_half(spec$EvoIso)
spec$n_individuals <- scale_half(log(spec$n_individuals))

Presence of disease models

Presence of disease base model

if (!file.exists("../fit_models/spec_disease_base.rds")) {

  spec_disease_m0 <- brm(Disease_YN ~ Host_order, 
    data=spec, family=bernoulli(link = "logit"), 
    prior=c(prior(normal(0, 1), class = Intercept),
                  prior(normal(0, 1.5), class = b)),
    iter=2000, thin=1,
    control=list(adapt_delta=0.8, max_treedepth=10), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_disease_m0, "../fit_models/spec_disease_base.rds")

} else { spec_disease_m0 <- readRDS("../fit_models/spec_disease_base.rds")}

model <- spec_disease_m0
# summary(model, prob=0.9)

# save memory
# rm(spec_disease_m0)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Mean CI CI_low CI_high Rhat
b_Intercept 0.07 0.90 -0.25 0.40 1.00
b_Host_orderChiroptera 1.02 0.90 0.53 1.52 1.00 
p1 <- cond_eff_plot(model)
ar1 <- intervals_plot(model)
plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

Sample size: 193

With no hierarchical effects and only the bat/rodent predictor, the model suggests that bats are more likely than rodents to suffer overt disease.

Presence of disease base model + phylogenetic effects

However, once include hierarchical effects for host and virus species, this effect goes away:

if (!file.exists("../fit_models/spec_disease_base_phylo.rds")) {

  spec_disease_m1 <- brm(Disease_YN ~ Host_order +
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=3000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.8, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_disease_m1, "../fit_models/spec_disease_base_phylo.rds")

} else { spec_disease_m1 <- readRDS("../fit_models/spec_disease_base_phylo.rds")}

model <- spec_disease_m1
# summary(model, prob=0.9)

# save memory
# rm(spec_disease_m1)

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.63 0.90 -0.68 1.94 1.00
b_Host_orderChiroptera fixed -0.06 0.90 -1.44 1.30 1.00
sd_Host_name__Intercept random 0.56 0.90 0.05 1.32 1.00 Host_name
sd_Host_Upham__Intercept random 0.81 0.90 0.08 1.78 1.00 Host_Upham
sd_Virus_ICTV__Intercept random 1.44 0.90 0.30 2.58 1.01 Virus_ICTV
sd_Virus_name__Intercept random 1.66 0.90 0.64 2.60 1.01 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-10,10)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

ggsave("../plots_tables/base_phylo_species_level_disease_comboplot.pdf", width=8, height=3)

Sample size: 193

Presence of disease full model

if (!file.exists("../fit_models/spec_disease_full.rds")) {

  spec_disease_full <- 
        brm(Disease_YN ~ Host_order + n_individuals + 
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.92, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_disease_full, "../fit_models/spec_disease_full.rds")

} else { spec_disease_full <- readRDS("../fit_models/spec_disease_full.rds")}

model <- spec_disease_full
# summary(model, prob=0.9)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.98 0.90 -0.35 2.29 1.00
b_Host_orderChiroptera fixed -9.55e-03 0.90 -1.54 1.50 1.00
b_n_individuals fixed 1.27 0.90 0.06 2.47 1.00
b_ses_pd fixed -0.69 0.90 -2.17 0.82 1.00
b_EvoIso fixed 0.39 0.90 -0.71 1.50 1.00
sd_Host_name__Intercept random 0.57 0.90 0.05 1.37 1.00 Host_name
sd_Host_Upham__Intercept random 0.80 0.90 0.07 1.82 1.00 Host_Upham
sd_PaperID__Intercept random 1.74 0.90 0.92 2.64 1.00 PaperID
sd_Virus_ICTV__Intercept random 1.17 0.90 0.12 2.40 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.14 0.90 0.15 2.23 1.00 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_disease_full_comboplot.pdf", width=8, height=3)

Sample size: 193

Mortality full model

if (!file.exists("../fit_models/spec_mort_full.rds")) {

  spec_mort_full <- 
        brm(Mortality_YN ~ Host_order + n_individuals + 
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_mort_full, "../fit_models/spec_mort_full.rds")

} else { spec_mort_full <- readRDS("../fit_models/spec_mort_full.rds")}

model <- spec_mort_full
# summary(model, prob=0.9)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.05 0.90 -1.13 1.20 1.00
b_Host_orderChiroptera fixed 0.49 0.90 -1.05 2.07 1.00
b_n_individuals fixed 2.34 0.90 1.00 3.71 1.00
b_ses_pd fixed -1.64 0.90 -2.92 -0.35 1.00
b_EvoIso fixed 1.03 0.90 -0.05 2.15 1.00
sd_Host_name__Intercept random 0.75 0.90 0.07 1.66 1.00 Host_name
sd_Host_Upham__Intercept random 0.93 0.90 0.12 1.93 1.00 Host_Upham
sd_PaperID__Intercept random 1.27 0.90 0.33 2.19 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.71 0.90 0.06 1.79 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.79 0.90 0.07 1.78 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of mortality") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_mort_full_comboplot.pdf", width=8, height=3)

Sample size: 154

Severity full model

if (!file.exists("../fit_models/spec_severity_full.rds")) {

  spec_severity_full <- 
        brm(severity_rank ~ Host_order + n_individuals + 
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=14), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_severity_full, "../fit_models/spec_severity_full.rds")

} else { spec_severity_full <- readRDS("../fit_models/spec_severity_full.rds")}

model <- spec_severity_full
# summary(model, prob=0.9)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -0.97 0.90 -1.88 -0.09 1.00
b_Intercept[2] fixed -0.53 0.90 -1.42 0.35 1.00
b_Intercept[3] fixed -0.25 0.90 -1.13 0.64 1.00
b_Intercept[4] fixed 0.49 0.90 -0.40 1.41 1.00
b_Host_orderChiroptera fixed -0.25 0.90 -1.49 0.97 1.00
b_n_individuals fixed 1.46 0.90 0.74 2.27 1.00
b_ses_pd fixed -0.66 0.90 -1.65 0.32 1.00
b_EvoIso fixed 0.32 0.90 -0.35 1.00 1.00
sd_Host_name__Intercept random 0.57 0.90 0.06 1.29 1.00 Host_name
sd_Host_Upham__Intercept random 0.80 0.90 0.09 1.64 1.00 Host_Upham
sd_PaperID__Intercept random 0.94 0.90 0.51 1.43 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.95 0.90 0.16 1.86 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.77 0.90 0.10 1.47 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of severity") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_severity_full_comboplot.pdf", width=8, height=3)

Sample size: 193

Joint results plots (species level models)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-12,12)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/species_level_conditionals_intervals_full.pdf", joint_plot, width=10.5, height=5.5)

Sample size: 193

Hierarchical effects plots

Species level disease full model

Species level mortality full model

Species level severity full model

Species level models (reservoir match)

Presence of disease (reservoir match)

if (!file.exists("../fit_models/spec_disease_resMatch.rds")) {

  spec_disease_resMatch <- brm(Disease_YN ~ Host_order + n_individuals + 
                      Reservoir_match + ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_disease_resMatch, "../fit_models/spec_disease_resMatch.rds")

} else { spec_disease_resMatch <- readRDS("../fit_models/spec_disease_resMatch.rds")}

model <- spec_disease_resMatch
# summary(model, prob=0.9)

# save memory
rm(spec_disease_resMatch)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.98 0.90 -0.59 2.56 1.00
b_Host_orderChiroptera fixed -0.41 0.90 -1.96 1.12 1.00
b_n_individuals fixed 1.02 0.90 -0.12 2.21 1.00
b_Reservoir_match fixed 0.09 0.90 -1.30 1.45 1.00
b_ses_pd fixed -0.27 0.90 -1.83 1.30 1.00
b_EvoIso fixed 0.26 0.90 -0.92 1.42 1.00
sd_Host_name__Intercept random 0.57 0.90 0.04 1.40 1.00 Host_name
sd_Host_Upham__Intercept random 0.80 0.90 0.07 1.82 1.00 Host_Upham
sd_PaperID__Intercept random 1.44 0.90 0.47 2.42 1.00 PaperID
sd_Virus_ICTV__Intercept random 1.50 0.90 0.21 2.82 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.22 0.90 0.16 2.36 1.00 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_disease_resMatch_comboplot.pdf", width=8, height=3)

Sample size: 176

Mortality (reservoir match)

if (!file.exists("../fit_models/spec_mort_resMatch.rds")) {

  spec_mort_resMatch <- brm(Mortality_YN ~ Host_order + n_individuals + 
                      Reservoir_match + ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_mort_resMatch, "../fit_models/spec_mort_resMatch.rds")

} else { spec_mort_resMatch <- readRDS("../fit_models/spec_mort_resMatch.rds")}

model <- spec_mort_resMatch
# summary(model, prob=0.90)

# save memory
rm(spec_mort_resMatch)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.12 0.90 -1.28 1.48 1.00
b_Host_orderChiroptera fixed 0.47 0.90 -1.10 2.03 1.00
b_n_individuals fixed 2.10 0.90 0.82 3.42 1.00
b_Reservoir_match fixed -0.20 0.90 -1.52 1.12 1.00
b_ses_pd fixed -1.38 0.90 -2.82 0.06 1.00
b_EvoIso fixed 0.94 0.90 -0.25 2.13 1.00
sd_Host_name__Intercept random 0.75 0.90 0.06 1.71 1.00 Host_name
sd_Host_Upham__Intercept random 0.93 0.90 0.10 1.94 1.00 Host_Upham
sd_PaperID__Intercept random 1.20 0.90 0.26 2.18 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.80 0.90 0.07 1.95 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.93 0.90 0.10 2.02 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of mortality") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_mortality_resMatch_comboplot.pdf", width=8, height=3)

Sample size: 139

Severity (reservoir match)

if (!file.exists("../fit_models/spec_severity_resMatch.rds")) {

  spec_severity_resMatch <- brm(severity_rank ~ Host_order + n_individuals + 
                      Reservoir_match + ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV),  
    data=spec, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_severity_resMatch, "../fit_models/spec_severity_resMatch.rds")

} else { spec_severity_resMatch <- readRDS("../fit_models/spec_severity_resMatch.rds")}

model <- spec_severity_resMatch
# summary(model, prob=0.90)

# save memory
rm(spec_severity_resMatch)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -1.02 0.90 -2.10 0.03 1.00
b_Intercept[2] fixed -0.67 0.90 -1.75 0.37 1.00
b_Intercept[3] fixed -0.36 0.90 -1.44 0.70 1.00
b_Intercept[4] fixed 0.39 0.90 -0.68 1.48 1.00
b_Host_orderChiroptera fixed -0.39 0.90 -1.63 0.87 1.00
b_n_individuals fixed 1.29 0.90 0.58 2.09 1.00
b_Reservoir_match fixed -0.10 0.90 -1.04 0.81 1.00
b_ses_pd fixed -0.35 0.90 -1.47 0.79 1.00
b_EvoIso fixed 0.20 0.90 -0.60 0.99 1.00
sd_Host_name__Intercept random 0.63 0.90 0.05 1.39 1.00 Host_name
sd_Host_Upham__Intercept random 0.76 0.90 0.09 1.59 1.00 Host_Upham
sd_PaperID__Intercept random 0.89 0.90 0.43 1.39 1.00 PaperID
sd_Virus_ICTV__Intercept random 1.15 0.90 0.19 2.15 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.90 0.90 0.15 1.72 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of severity") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_severity_resMatch_comboplot.pdf", width=8, height=3)

Sample size: 176

Joint results plots (species level reservoir match models)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-12,12)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/species_level_conditionals_intervals_resMatch.pdf", joint_plot, width=10.5, height=5.5)

Sample size: 176

Species level models (excluding Lyssaviruses)

Presence of disease (excluding Lyssaviruses)

lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
spec_nolys <- spec[!spec$Virus_ICTV%in%lyssaviruses,]

if (!file.exists("../fit_models/spec_disease_nolys.rds")) {

  spec_disease_nolys <- brm(Disease_YN ~ Host_order + n_individuals + 
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec_nolys, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_disease_nolys, "../fit_models/spec_disease_nolys.rds")

} else { spec_disease_nolys <- readRDS("../fit_models/spec_disease_nolys.rds")}

model <- spec_disease_nolys
# summary(model, prob=0.9)

# # save memory
rm(spec_disease_nolys)

# # pp <- brms::posterior_predict(model)
# # ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

# param_tab(model)

p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditonal log-odds of disease") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_disease_nolys.pdf", width=8, height=3)

Sample size: 148

Mortality (excluding Lyssaviruses)

if (!file.exists("../fit_models/spec_mort_nolys.rds")) {

  spec_mort_nolys <- brm(Mortality_YN ~ Host_order + n_individuals +  
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec_nolys, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_mort_nolys, "../fit_models/spec_mort_nolys.rds")

} else { spec_mort_nolys <- readRDS("../fit_models/spec_mort_nolys.rds")}

model <- spec_mort_nolys
# summary(model, prob=0.9)

# save memory
rm(spec_mort_nolys)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.14 0.90 -0.97 1.26 1.00
b_Host_orderChiroptera fixed -0.10 0.90 -1.69 1.53 1.00
b_n_individuals fixed 1.94 0.90 0.49 3.43 1.00
b_ses_pd fixed -1.62 0.90 -3.18 -0.09 1.00
b_EvoIso fixed 1.33 0.90 0.16 2.55 1.00
sd_Host_name__Intercept random 0.66 0.90 0.05 1.55 1.00 Host_name
sd_Host_Upham__Intercept random 0.75 0.90 0.07 1.73 1.00 Host_Upham
sd_PaperID__Intercept random 1.39 0.90 0.46 2.34 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.67 0.90 0.05 1.71 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.76 0.90 0.06 1.76 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of mortality") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_mortality_nolys_comboplot.pdf", width=8, height=3)

Sample size: 113

Severity (excluding Lyssaviruses)

if (!file.exists("../fit_models/spec_severity_nolys.rds")) {

  spec_severity_nolys <- brm(severity_rank ~ Host_order + n_individuals +
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV),  
    data=spec_nolys, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_severity_nolys, "../fit_models/spec_severity_nolys.rds")

} else { spec_severity_nolys <- readRDS("../fit_models/spec_severity_nolys.rds")}

model <- spec_severity_nolys
# summary(model, prob=0.9)

# save memory
rm(spec_severity_nolys)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -1.06 0.90 -1.91 -0.23 1.00
b_Intercept[2] fixed -0.58 0.90 -1.40 0.24 1.00
b_Intercept[3] fixed -0.29 0.90 -1.10 0.52 1.00
b_Intercept[4] fixed 0.35 0.90 -0.46 1.19 1.00
b_Host_orderChiroptera fixed -0.41 0.90 -1.66 0.93 1.00
b_n_individuals fixed 1.23 0.90 0.41 2.07 1.00
b_ses_pd fixed -1.04 0.90 -2.15 0.02 1.00
b_EvoIso fixed 0.46 0.90 -0.20 1.14 1.00
sd_Host_name__Intercept random 0.56 0.90 0.05 1.34 1.01 Host_name
sd_Host_Upham__Intercept random 0.85 0.90 0.10 1.75 1.00 Host_Upham
sd_PaperID__Intercept random 0.94 0.90 0.49 1.44 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.45 0.90 0.04 1.16 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.84 0.90 0.16 1.53 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of severity") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_severity_nolys_comboplot.pdf", width=8, height=3)

Sample size: 148

Joint results plots (excluding Lyssaviruses)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-12,12)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/species_level_conditionals_intervals_nolys.pdf", joint_plot, width=10.5, height=5.5)

Species level models (only heterologous hosts)

Presence of disease (heterologous hosts)

spec_heteros <- spec[spec$Reservoir_match%in%0,]

if (!file.exists("../fit_models/spec_disease_heteroHosts.rds")) {

  spec_disease_heteros <- brm(Disease_YN ~ Host_order + n_individuals + 
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec_heteros, family=bernoulli(link = "logit"),
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=10), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_disease_heteros, "../fit_models/spec_disease_heteroHosts.rds")

} else { spec_disease_heteros <- readRDS("../fit_models/spec_disease_heteroHosts.rds")}

model <- spec_disease_heteros
# summary(model, prob=0.9)

# # save memory
rm(spec_disease_heteros)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

# param_tab(model)

p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditonal log-odds of disease") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(ar1)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_disease_heteroHosts_comboplot.pdf", width=8, height=3)

Sample size: 81

Mortality (heterologous hosts)

if (!file.exists("../fit_models/spec_mort_heteroHosts.rds")) {

  spec_mort_heteros <- brm(Mortality_YN ~ Host_order + n_individuals +  
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV), 
    data=spec_heteros, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=10), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_mort_heteros, "../fit_models/spec_mort_heteroHosts.rds")

} else { spec_mort_heteros <- readRDS("../fit_models/spec_mort_heteroHosts.rds")}

model <- spec_mort_heteros
# summary(model, prob=0.9)

# save memory
rm(spec_mort_heteros)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.40 0.90 -0.87 1.66 1.00
b_Host_orderChiroptera fixed 0.16 0.90 -1.51 1.84 1.00
b_n_individuals fixed 2.05 0.90 0.35 3.74 1.00
b_ses_pd fixed -0.91 0.90 -2.78 0.94 1.00
b_EvoIso fixed 0.23 0.90 -1.28 1.74 1.00
sd_Host_name__Intercept random 0.92 0.90 0.08 2.12 1.00 Host_name
sd_Host_Upham__Intercept random 0.69 0.90 0.06 1.70 1.00 Host_Upham
sd_PaperID__Intercept random 0.90 0.90 0.08 2.06 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.64 0.90 0.05 1.64 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.91 0.90 0.08 2.03 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of mortality") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(ar2)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_mortality_heteroHosts_comboplot.pdf", width=8, height=3)

Sample size: 61

Severity (heterologous hosts)

if (!file.exists("../fit_models/spec_severity_heteroHosts.rds")) {

  spec_severity_heteros <- brm(severity_rank ~ Host_order + n_individuals +
                      ses_pd + EvoIso + 
                      (1|PaperID) + 
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV),  
    data=spec_heteros, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=10), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(spec_severity_heteros, "../fit_models/spec_severity_heteroHosts.rds")

} else { spec_severity_heteros <- readRDS("../fit_models/spec_severity_heteroHosts.rds")}

# pairs(spec_severity_heteros)
model <- spec_severity_heteros
# summary(model, prob=0.9)

# save memory
rm(spec_severity_heteros)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -1.15 0.90 -2.18 -0.15 1.00
b_Intercept[2] fixed -1.00 0.90 -2.03 -7.83e-03 1.00
b_Intercept[3] fixed -0.57 0.90 -1.57 0.40 1.00
b_Intercept[4] fixed 0.23 0.90 -0.77 1.24 1.00
b_Host_orderChiroptera fixed -0.40 0.90 -1.74 0.97 1.00
b_n_individuals fixed 1.67 0.90 0.44 3.01 1.00
b_ses_pd fixed -1.10 0.90 -2.77 0.54 1.00
b_EvoIso fixed 0.05 0.90 -1.20 1.35 1.00
sd_Host_name__Intercept random 1.42 0.90 0.47 2.39 1.00 Host_name
sd_Host_Upham__Intercept random 0.55 0.90 0.04 1.40 1.00 Host_Upham
sd_PaperID__Intercept random 0.58 0.90 0.05 1.35 1.00 PaperID
sd_Virus_ICTV__Intercept random 0.70 0.90 0.05 1.75 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.34 0.90 0.43 2.30 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + ylab("Conditional log-odds of severity") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/species_level_severity_heteroHosts_comboplot.pdf", width=8, height=3)

Sample size: 81

Joint results plots (heterologous hosts)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-12,12)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/species_level_conditionals_intervals_heteroHosts.pdf", joint_plot, width=10.5, height=5.5)

Sample size: 81

Individual level models

Modelling dose differences between bats and rodents

Parameter Effects Component Mean CI CI_low CI_high Rhat Group
b_Intercept fixed conditional 1.22 0.90 -0.64 3.09 1.00
b_Host_orderChiroptera fixed conditional 0.64 0.90 -1.30 2.43 1.00
sd_Dose_unit__Intercept random conditional 1.98 0.90 1.11 3.02 1.00 Dose_unit
sd_Host_name__Intercept random conditional 2.05 0.90 1.52 2.58 1.00 Host_name
sd_Host_Upham__Intercept random conditional 1.87 0.90 0.74 2.95 1.00 Host_Upham
sd_Route_type__Intercept random conditional 1.56 0.90 0.82 2.58 1.00 Route_type
sd_Virus_ICTV__Intercept random conditional 1.05 0.90 0.13 2.22 1.00 Virus_ICTV
sd_Virus_name__Intercept random conditional 2.08 0.90 1.61 2.60 1.00 Virus_name
sigma fixed sigma 1.82 0.90 1.76 1.88 1.00

Presence of disease

Presence of disease - base model

Base model with no hierarchical effects:

if (!file.exists("../fit_models/indiv_disease_base.rds")) {

  disease_m0 <- brm(Disease_YN ~ Host_order, 
    data=dat_dis, family=bernoulli(link = "logit"), 
    prior=c(prior(normal(0, 1), class = Intercept),
                  prior(normal(0, 1.5), class = b)),
    iter=4000, thin=1,
    control=list(adapt_delta=0.8, max_treedepth=10), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(disease_m0, "../fit_models/indiv_disease_base.rds")

} else { disease_m0 <- readRDS("../fit_models/indiv_disease_base.rds")}

model <- disease_m0
# summary(model, prob=0.9)

# save memory
# rm(disease_m0)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Mean CI CI_low CI_high Rhat
b_Intercept -0.57 0.90 -0.71 -0.44 1.00
b_Host_orderChiroptera 0.40 0.90 0.24 0.56 1.00 
p1 <- cond_eff_plot(model) + ylab("Conditional log-odds of disease") + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

Sample size: 2444

With no hierarchical effects and only the bat / rodent binary predictor, the model suggests that bats are more likely than rodents to suffer overt disease.

However, once include hierarchical effects for host and virus species, this effect goes away:

Presence of disease - base model + phylogenetic effects

if (!file.exists("../fit_models/indiv_disease_base_phylo.rds")) {

  disease_m1 <- brm(Disease_YN ~ Host_order +
                    (1|Host_name) + (1|Virus_name) +
                    (1|Host_Upham) + (1|Virus_ICTV), 
    data=dat_dis, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.99, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(disease_m1, "../fit_models/indiv_disease_base_phylo.rds")

} else { disease_m1 <- readRDS("../fit_models/indiv_disease_base_phylo.rds")}

model <- disease_m1
# summary(model, prob=0.9)

# save memory
# rm(disease_m1)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.47 0.90 -1.24 2.18 1.00
b_Host_orderChiroptera fixed -0.69 0.90 -2.37 1.04 1.00
sd_Host_name__Intercept random 3.41 0.90 2.69 4.18 1.00 Host_name
sd_Host_Upham__Intercept random 1.15 0.90 0.11 2.50 1.00 Host_Upham
sd_Virus_ICTV__Intercept random 1.95 0.90 0.76 3.12 1.00 Virus_ICTV
sd_Virus_name__Intercept random 2.58 0.90 1.80 3.41 1.00 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))
ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_disease_base_phylo_comboplot.pdf", width=8, height=3)

Sample size: 2444

Presence of disease - Full model

Full model adjusts for evolutionary predictors, dose, and inoculation route

if (!file.exists("../fit_models/indiv_disease_full.rds")) {

  disease_full <- brm(Disease_YN ~ Host_order + 
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_dis, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(disease_full, "../fit_models/indiv_disease_full.rds")

} else { disease_full <- readRDS("../fit_models/indiv_disease_full.rds")}

model <- disease_full
# summary(model, prob=0.9)

# save memory
# rm(disease_full)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.18 0.90 -1.64 2.01 1.00
b_Host_orderChiroptera fixed -0.54 0.90 -2.20 1.14 1.00
b_ses_pd fixed -1.46 0.90 -2.97 0.04 1.00
b_EvoIso fixed 0.70 0.90 -0.55 1.99 1.00
b_Dose_mass fixed 0.71 0.90 0.26 1.18 1.00
sd_Dose_unit__Intercept random 2.30 0.90 1.48 3.24 1.00 Dose_unit
sd_Host_name__Intercept random 2.77 0.90 2.10 3.47 1.00 Host_name
sd_Host_Upham__Intercept random 1.01 0.90 0.08 2.36 1.00 Host_Upham
sd_Route_type__Intercept random 1.62 0.90 0.98 2.45 1.00 Route_type
sd_Virus_ICTV__Intercept random 3.17 0.90 2.06 4.21 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.19 0.90 0.14 2.39 1.00 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank()) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_disease_full_comboplot.pdf", width=8, height=3)

Sample size: 1500

Mortality full model

dat_mort <- dat
dat_mort <- dat_mort[dat_mort$Susceptible_YN==1,]
dat_mort <- dat_mort[!is.na(dat_mort$Susceptible_YN),]

dat_mort <- dat_mort[!is.na(dat_mort$Mortality_YN),]

dat_mort$ses_pd <- scale_half(dat_mort$ses_pd)
dat_mort$EvoIso <- scale_half(dat_mort$EvoIso)

dat_mort <- dat_mort %>% group_by(Dose_unit) %>% mutate(Dose_mass = scale_half(log(Dose_mass)))

dat_mort$Host_name <- dat_mort$Host_Upham
dat_mort$Virus_name <- dat_mort$Virus_ICTV
if (!file.exists("../fit_models/indiv_mort_full.rds")) {

  mort_full <- brm(Mortality_YN ~ Host_order +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_mort, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, 
    cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(mort_full, "../fit_models/indiv_mort_full.rds")

} else { mort_full <- readRDS("../fit_models/indiv_mort_full.rds")}

model <- mort_full
# summary(model, prob=0.9)

# save memory
# rm(mort_full)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed -1.59 0.90 -3.22 0.10 1.00
b_Host_orderChiroptera fixed 1.31 0.90 -0.40 3.01 1.00
b_ses_pd fixed -2.25 0.90 -3.32 -1.18 1.00
b_EvoIso fixed 0.43 0.90 -0.79 1.74 1.00
b_Dose_mass fixed 5.75e-03 0.90 -0.52 0.54 1.00
sd_Dose_unit__Intercept random 1.75 0.90 1.05 2.65 1.00 Dose_unit
sd_Host_name__Intercept random 2.28 0.90 1.66 2.92 1.00 Host_name
sd_Host_Upham__Intercept random 0.88 0.90 0.08 2.03 1.00 Host_Upham
sd_Route_type__Intercept random 1.92 0.90 1.15 2.87 1.00 Route_type
sd_Virus_ICTV__Intercept random 0.59 0.90 0.04 1.50 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.59 0.90 0.05 1.42 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,18)) + 
      ylab("Conditional log-odds of mortality") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_mortality_full_comboplot.pdf", width=8, height=3)

Sample size: 927

Severity full model

Fitting an ordinal regression (a.k.a. cumulative link) model for disease severity

dat_severity <- dat
dat_severity <- dat_severity[dat_severity$Susceptible_YN==1,]
dat_severity <- dat_severity[!is.na(dat_severity$Susceptible_YN),]

dat_severity <- dat_severity[!is.na(dat_severity$severity_rank),]

dat_severity$ses_pd <- scale_half(dat_severity$ses_pd)
dat_severity$EvoIso <- scale_half(dat_severity$EvoIso)

dat_severity <- dat_severity %>% group_by(Dose_unit) %>% mutate(Dose_mass = scale_half(log(Dose_mass)))

dat_severity$Host_name <- dat_severity$Host_Upham
dat_severity$Virus_name <- dat_severity$Virus_ICTV
if (!file.exists("../fit_models/indiv_severity_full.rds")) {

  severity_full <- brm(severity_rank ~ Host_order + 
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_severity, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.98, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(severity_full, "../fit_models/indiv_severity_full.rds")

} else { severity_full <- readRDS("../fit_models/indiv_severity_full.rds")}

model <- severity_full
# summary(model, prob=0.9)

# save memory
# rm(severity_full)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -0.52 0.90 -1.64 0.64 1.00
b_Intercept[2] fixed 0.02 0.90 -1.10 1.19 1.00
b_Intercept[3] fixed 0.32 0.90 -0.80 1.48 1.00
b_Intercept[4] fixed 0.75 0.90 -0.37 1.92 1.00
b_Host_orderChiroptera fixed -0.06 0.90 -1.32 1.30 1.00
b_ses_pd fixed -1.16 0.90 -2.05 -0.27 1.00
b_EvoIso fixed 0.22 0.90 -0.43 0.91 1.00
b_Dose_mass fixed 0.60 0.90 0.40 0.80 1.00
sd_Dose_unit__Intercept random 1.72 0.90 1.04 2.61 1.00 Dose_unit
sd_Host_name__Intercept random 1.48 0.90 0.99 1.95 1.00 Host_name
sd_Host_Upham__Intercept random 0.90 0.90 0.08 1.97 1.01 Host_Upham
sd_Route_type__Intercept random 1.00 0.90 0.51 1.75 1.00 Route_type
sd_Virus_ICTV__Intercept random 2.01 0.90 1.39 2.67 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.40 0.90 0.03 0.98 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-12,12)) + 
      ylab("Conditional log-odds of severity") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_severity_full_comboplot.pdf", width=8, height=3)

Sample size: 1500

Joint results plots (individual level models)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-15,15)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/individual_level_conditionals_intervals_full.pdf", joint_plot, width=10.5, height=5.5)

Hierarchical effects plots

Individual level presence of disease full model

Individual level mortality full model

Individual level severity full model

Individual level models (reservoir match)

Presence of disease (reservoir match)

if (!file.exists("../fit_models/indiv_disease_resMatch.rds")) {

  disease_resMatch <- brm(Disease_YN ~ Host_order + Reservoir_match +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_dis, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(disease_resMatch, "../fit_models/indiv_disease_resMatch.rds")

} else { disease_resMatch <- readRDS("../fit_models/indiv_disease_resMatch.rds")}

model <- disease_resMatch
# summary(model, prob=0.9)

# save memory
# rm(disease_resMatch)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed -0.50 0.90 -2.68 1.67 1.00
b_Host_orderChiroptera fixed -0.54 0.90 -2.19 1.13 1.00
b_Reservoir_match fixed 1.09 0.90 -0.35 2.57 1.00
b_ses_pd fixed -1.26 0.90 -2.85 0.28 1.00
b_EvoIso fixed 0.58 0.90 -0.77 1.87 1.00
b_Dose_mass fixed 0.67 0.90 0.21 1.12 1.00
sd_Dose_unit__Intercept random 2.18 0.90 1.38 3.11 1.00 Dose_unit
sd_Host_name__Intercept random 2.60 0.90 1.93 3.33 1.00 Host_name
sd_Host_Upham__Intercept random 0.90 0.90 0.07 2.19 1.00 Host_Upham
sd_Route_type__Intercept random 1.67 0.90 1.01 2.50 1.00 Route_type
sd_Virus_ICTV__Intercept random 2.59 0.90 0.87 3.97 1.01 Virus_ICTV
sd_Virus_name__Intercept random 1.49 0.90 0.20 2.76 1.01 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_disease_resMatch_comboplot.pdf", width=8, height=3)

Sample size: 1391

Mortality (reservoir match)

if (!file.exists("../fit_models/indiv_mortality_resMatch.rds")) {

  mortality_resMatch <- brm(Mortality_YN ~ Host_order + Reservoir_match +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_mort, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(mortality_resMatch, "../fit_models/indiv_mortality_resMatch.rds")

} else { mortality_resMatch <- readRDS("../fit_models/indiv_mortality_resMatch.rds")}

model <- mortality_resMatch
# summary(model, prob=0.9)

# save memory
# rm(mortality_resMatch)

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed -1.62 0.90 -3.45 0.20 1.00
b_Host_orderChiroptera fixed 1.03 0.90 -0.59 2.64 1.00
b_Reservoir_match fixed 0.29 0.90 -1.06 1.69 1.00
b_ses_pd fixed -1.97 0.90 -3.15 -0.71 1.00
b_EvoIso fixed 0.51 0.90 -0.80 1.90 1.00
b_Dose_mass fixed 0.02 0.90 -0.52 0.57 1.00
sd_Dose_unit__Intercept random 1.70 0.90 1.00 2.56 1.00 Dose_unit
sd_Host_name__Intercept random 2.19 0.90 1.60 2.82 1.00 Host_name
sd_Host_Upham__Intercept random 0.75 0.90 0.06 1.81 1.00 Host_Upham
sd_Route_type__Intercept random 1.97 0.90 1.18 2.95 1.00 Route_type
sd_Virus_ICTV__Intercept random 0.66 0.90 0.05 1.66 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.74 0.90 0.07 1.69 1.00 Virus_name 
# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Mortality_YN, pp[1:500, ])

p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) + 
      ylab("Conditional log-odds of mortality") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_mortality_resMatch_comboplot.pdf", width=8, height=3)

Sample size: 851

Severity (reservoir match)

if (!file.exists("../fit_models/indiv_severity_resMatch.rds")) {

  severity_resMatch <- brm(severity_rank ~ Host_order + Reservoir_match +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_severity, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.98, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(severity_resMatch, "../fit_models/indiv_severity_resMatch.rds")

} else { severity_resMatch <- readRDS("../fit_models/indiv_severity_resMatch.rds")}

model <- severity_resMatch
# summary(model, prob=0.9)

# save memory
# rm(severity_resMatch)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$severity_rank, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed 0.01 0.90 -1.37 1.40 1.00
b_Intercept[2] fixed 0.50 0.90 -0.89 1.89 1.00
b_Intercept[3] fixed 0.80 0.90 -0.59 2.20 1.00
b_Intercept[4] fixed 1.24 0.90 -0.14 2.63 1.00
b_Host_orderChiroptera fixed -8.57e-03 0.90 -1.22 1.27 1.00
b_Reservoir_match fixed 0.72 0.90 -0.19 1.66 1.00
b_ses_pd fixed -1.02 0.90 -2.08 -0.04 1.00
b_EvoIso fixed 0.16 0.90 -0.64 0.98 1.00
b_Dose_mass fixed 0.59 0.90 0.39 0.79 1.00
sd_Dose_unit__Intercept random 1.69 0.90 1.01 2.56 1.00 Dose_unit
sd_Host_name__Intercept random 1.43 0.90 0.98 1.88 1.00 Host_name
sd_Host_Upham__Intercept random 0.72 0.90 0.05 1.81 1.00 Host_Upham
sd_Route_type__Intercept random 1.00 0.90 0.51 1.76 1.00 Route_type
sd_Virus_ICTV__Intercept random 1.78 0.90 0.68 2.71 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.74 0.90 0.09 1.54 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,15)) + 
      ylab("Conditional log-odds of severity") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_severity_resMatch_comboplot.pdf", width=8, height=3)

Sample size: 1391

Joint results plots (reservoir match)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-18,18)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/individual_level_conditionals_intervals_resMatch.pdf", joint_plot, width=10.5, height=5.5)

Individual level models (excluding Lyssaviruses)

Presence of disease (excluding Lyssaviruses)

lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
dat_dis_nolys <- dat_dis[!dat_dis$Virus_ICTV%in%lyssaviruses,]

if (!file.exists("../fit_models/indiv_disease_nolys.rds")) {

  disease_nolys <- brm(Disease_YN ~ Host_order + 
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_dis_nolys, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.9, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(disease_nolys, "../fit_models/indiv_disease_nolys.rds")

} else { disease_nolys <- readRDS("../fit_models/indiv_disease_nolys.rds")}

model <- disease_nolys
# summary(model, prob=0.9)

# save memory
# rm(disease_nolys)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed 0.23 0.90 -1.61 2.06 1.00
b_Host_orderChiroptera fixed -0.62 0.90 -2.34 1.20 1.00
b_ses_pd fixed -1.15 0.90 -3.16 0.82 1.00
b_EvoIso fixed 0.77 0.90 -0.73 2.33 1.00
b_Dose_mass fixed 1.05 0.90 0.37 1.76 1.00
sd_Dose_unit__Intercept random 1.46 0.90 0.31 2.65 1.00 Dose_unit
sd_Host_name__Intercept random 3.26 0.90 2.53 4.05 1.00 Host_name
sd_Host_Upham__Intercept random 1.04 0.90 0.08 2.49 1.00 Host_Upham
sd_Route_type__Intercept random 1.24 0.90 0.50 2.18 1.00 Route_type
sd_Virus_ICTV__Intercept random 1.77 0.90 0.21 3.43 1.00 Virus_ICTV
sd_Virus_name__Intercept random 2.87 0.90 1.69 3.97 1.00 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,18)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_disease_nolys_comboplot.pdf", width=8, height=3)

Sample size: 1425

Mortality (excluding Lyssaviruses)

lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
dat_mort_nolys <- dat_mort[!dat_mort$Virus_ICTV%in%lyssaviruses,]

if (!file.exists("../fit_models/indiv_mort_nolys.rds")) {

  mort_nolys <- brm(Mortality_YN ~ Host_order + 
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_mort_nolys, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(mort_nolys, "../fit_models/indiv_mort_nolys.rds")

} else { mort_nolys <- readRDS("../fit_models/indiv_mort_nolys.rds")}

model <- mort_nolys
# summary(model, prob=0.9)

# save memory
# rm(mort_nolys)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed -0.91 0.90 -2.56 0.78 1.00
b_Host_orderChiroptera fixed 0.19 0.90 -1.67 2.09 1.00
b_ses_pd fixed -1.55 0.90 -3.16 0.13 1.00
b_EvoIso fixed 1.40 0.90 -0.25 3.11 1.00
b_Dose_mass fixed 0.28 0.90 -0.60 1.18 1.00
sd_Dose_unit__Intercept random 1.32 0.90 0.18 2.55 1.00 Dose_unit
sd_Host_name__Intercept random 2.18 0.90 1.45 2.95 1.00 Host_name
sd_Host_Upham__Intercept random 1.05 0.90 0.09 2.37 1.00 Host_Upham
sd_Route_type__Intercept random 1.56 0.90 0.65 2.67 1.00 Route_type
sd_Virus_ICTV__Intercept random 0.84 0.90 0.07 2.00 1.00 Virus_ICTV
sd_Virus_name__Intercept random 0.96 0.90 0.10 2.02 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,10)) + 
      ylab("Conditional log-odds of mortality") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_mortality_nolys_comboplot.pdf", width=8, height=3)

Sample size: 742

Severity (excluding Lyssaviruses)

lyssaviruses <- viruses$Virus_ICTV[viruses$genus=="Lyssavirus"] %>% unique()
dat_severity_nolys <- dat_severity[!dat_severity$Virus_ICTV%in%lyssaviruses,]

if (!file.exists("../fit_models/indiv_severity_nolys.rds")) {

  severity_nolys <- brm(severity_rank ~ Host_order + 
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_severity_nolys, family=cumulative(probit),  
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(severity_nolys, "../fit_models/indiv_severity_nolys.rds")

} else { severity_nolys <- readRDS("../fit_models/indiv_severity_nolys.rds")}

model <- severity_nolys
# summary(model, prob=0.9)

# save memory
# rm(severity_nolys)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -1.17 0.90 -2.45 0.11 1.00
b_Intercept[2] fixed -0.34 0.90 -1.62 0.94 1.00
b_Intercept[3] fixed 0.08 0.90 -1.20 1.37 1.00
b_Intercept[4] fixed 0.65 0.90 -0.63 1.93 1.00
b_Host_orderChiroptera fixed -0.28 0.90 -1.75 1.32 1.00
b_ses_pd fixed -1.23 0.90 -2.69 0.25 1.00
b_EvoIso fixed 0.50 0.90 -0.42 1.41 1.00
b_Dose_mass fixed 0.42 0.90 0.16 0.69 1.00
sd_Dose_unit__Intercept random 1.28 0.90 0.38 2.24 1.00 Dose_unit
sd_Host_name__Intercept random 1.98 0.90 1.38 2.58 1.00 Host_name
sd_Host_Upham__Intercept random 1.08 0.90 0.09 2.42 1.01 Host_Upham
sd_Route_type__Intercept random 1.35 0.90 0.65 2.23 1.00 Route_type
sd_Virus_ICTV__Intercept random 1.74 0.90 0.43 2.82 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.11 0.90 0.14 2.10 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,10)) + 
      ylab("Conditional log-odds of severity") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_severity_nolys_comboplot.pdf", width=8, height=3)

Sample size: 1425

Joint results plots (excluding Lyssaviruses)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-18,10)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/individual_level_conditionals_intervals_nolys.pdf", joint_plot, width=10.5, height=5.5)

Individual level models (only heterologous hosts)

Presence of disease (heterologous hosts)

dat_dis_heteros <- dat_dis[dat_dis$Reservoir_match%in%0,]

if (!file.exists("../fit_models/indiv_disease_heteroHosts.rds")) {

  disease_heteros <- brm(Disease_YN ~ Host_order +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_dis_heteros, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(disease_heteros, "../fit_models/indiv_disease_heteroHosts.rds")

} else { disease_heteros <- readRDS("../fit_models/indiv_disease_heteroHosts.rds")}

model <- disease_heteros
# summary(model, prob=0.9)

# save memory
# rm(disease_heteros)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed -0.17 0.90 -2.12 1.78 1.00
b_Host_orderChiroptera fixed -0.67 0.90 -2.42 1.08 1.00
b_ses_pd fixed -0.62 0.90 -2.78 1.51 1.00
b_EvoIso fixed 1.16 0.90 -0.39 2.89 1.00
b_Dose_mass fixed 1.02 0.90 0.30 1.78 1.00
sd_Dose_unit__Intercept random 1.17 0.90 0.17 2.31 1.00 Dose_unit
sd_Host_name__Intercept random 2.66 0.90 1.94 3.45 1.00 Host_name
sd_Host_Upham__Intercept random 0.78 0.90 0.06 1.93 1.00 Host_Upham
sd_Route_type__Intercept random 1.60 0.90 0.84 2.55 1.00 Route_type
sd_Virus_ICTV__Intercept random 1.05 0.90 0.11 2.26 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.25 0.90 0.16 2.44 1.00 Virus_name 
p1 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-15,12)) + 
      ylab("Conditional log-odds of disease") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10)) + 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar1 <- intervals_plot(model)

plot_grid(p1, ar1, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_disease_heteroHosts_comboplot.pdf", width=8, height=3)

Sample size: 872

Mortality (heterologous hosts)

dat_mort_heteros <- dat_mort[dat_mort$Reservoir_match%in%0,]

if (!file.exists("../fit_models/indiv_mort_heteroHosts.rds")) {

  mort_heteros <- brm(Mortality_YN ~ Host_order  +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_mort_heteros, family=bernoulli(link = "logit"), 
    prior=custom_prior,
    iter=5000, thin=1, 
    cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.95, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(mort_heteros, "../fit_models/indiv_mort_heteroHosts.rds")

} else { mort_heteros <- readRDS("../fit_models/indiv_mort_heteroHosts.rds")}

model <- mort_heteros
# summary(model, prob=0.9)

# save memory
# rm(mort_heteros)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept fixed -1.06 0.90 -2.83 0.72 1.00
b_Host_orderChiroptera fixed 0.11 0.90 -1.71 1.91 1.00
b_ses_pd fixed 0.18 0.90 -1.98 2.30 1.00
b_EvoIso fixed 1.06 0.90 -0.79 3.07 1.00
b_Dose_mass fixed 0.37 0.90 -0.61 1.36 1.00
sd_Dose_unit__Intercept random 1.05 0.90 0.10 2.34 1.00 Dose_unit
sd_Host_name__Intercept random 1.45 0.90 0.51 2.40 1.00 Host_name
sd_Host_Upham__Intercept random 0.92 0.90 0.09 2.10 1.00 Host_Upham
sd_Route_type__Intercept random 1.78 0.90 0.93 2.80 1.00 Route_type
sd_Virus_ICTV__Intercept random 0.85 0.90 0.07 2.05 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.45 0.90 0.31 2.57 1.00 Virus_name 
p2 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-18,12)) + 
      ylab("Conditional log-odds of mortality") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar2 <- intervals_plot(model)

plot_grid(p2, ar2, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_mortality_heteroHosts_comboplot.pdf", width=8, height=3)

Sample size: 358

Severity (heterologous hosts)

dat_severity_heteros <- dat_severity[dat_severity$Reservoir_match%in%0,]

if (!file.exists("../fit_models/indiv_severity_heteroHosts.rds")) {

  severity_heteros <- brm(severity_rank ~ Host_order  +
                      ses_pd + EvoIso +  
                       Dose_mass + (1|Dose_unit) +    
                      (1|Host_name) + (1|Virus_name) +
                      (1|Host_Upham) + (1|Virus_ICTV) + 
                      (1|Route_type), 
    data=dat_severity_heteros, family=cumulative(probit), 
    prior=custom_prior,
    iter=5000, thin=1, cov_ranef = list(Host_Upham = host_cov, Virus_ICTV=virus_cov),
    control=list(adapt_delta=0.98, max_treedepth=12), cores=4, 
    save_pars=save_pars(all=TRUE))

  saveRDS(severity_heteros, "../fit_models/indiv_severity_heteroHosts.rds")

} else { severity_heteros <- readRDS("../fit_models/indiv_severity_heteroHosts.rds")}

model <- severity_heteros
# summary(model, prob=0.9)

# save memory
# rm(severity_heteros)

# pp <- brms::posterior_predict(model)
# ppc_dens_overlay(y=model$data$Disease_YN, pp[1:500, ])

param_tab(model)
Parameter Effects Mean CI CI_low CI_high Rhat Group
b_Intercept[1] fixed -0.47 0.90 -1.80 0.87 1.00
b_Intercept[2] fixed -0.27 0.90 -1.61 1.06 1.00
b_Intercept[3] fixed 0.08 0.90 -1.26 1.41 1.00
b_Intercept[4] fixed 0.65 0.90 -0.68 1.98 1.00
b_Host_orderChiroptera fixed -0.39 0.90 -1.73 0.98 1.00
b_ses_pd fixed -0.60 0.90 -2.45 1.21 1.00
b_EvoIso fixed 0.61 0.90 -0.64 2.02 1.00
b_Dose_mass fixed 0.47 0.90 0.16 0.78 1.00
sd_Dose_unit__Intercept random 0.94 0.90 0.13 1.89 1.00 Dose_unit
sd_Host_name__Intercept random 1.30 0.90 0.75 1.85 1.00 Host_name
sd_Host_Upham__Intercept random 0.75 0.90 0.07 1.74 1.00 Host_Upham
sd_Route_type__Intercept random 1.57 0.90 0.88 2.46 1.00 Route_type
sd_Virus_ICTV__Intercept random 0.93 0.90 0.10 2.03 1.00 Virus_ICTV
sd_Virus_name__Intercept random 1.17 0.90 0.31 2.04 1.00 Virus_name 
p3 <- cond_eff_plot(model) + scale_y_continuous(lim=c(-10,5)) + 
      ylab("Conditional log-odds of severity") + xlab("") + 
      theme(axis.title=element_text(size=9),axis.text.x=element_text(size=10), 
            panel.grid.major = element_blank(), 
            panel.grid.minor = element_blank())+ 
      scale_x_discrete(limits = c("Chiroptera","Rodentia"))

ar3 <- intervals_plot(model)

plot_grid(p3, ar3, align="h", rel_widths =c(0.7,1))

# ggsave("../plots_tables/individual_level_severity_heteroHosts_comboplot.pdf", width=8, height=3)

Sample size: 872

Joint results plots (individual level heterologous host models)

## Joint plot
global_ylim <-  scale_y_continuous(lim=c(-15,10)) 
padding <- theme(plot.margin = unit(c(t=0, r=-3, b=0, l=0), "cm"))

joint_plot <- joint_conditionals_intervals(p1,p2,p3,global_ylim,padding,ar1,ar2,ar3)

ggsave("../plots_tables/individual_level_conditionals_intervals_heteroHosts.pdf", joint_plot, width=10.5, height=5.5)

Joint results plots - species and individual models combined

s_dis <- readRDS("../fit_models/spec_disease_full.rds")
s_mort <- readRDS("../fit_models/spec_mort_full.rds")
s_sev <- readRDS("../fit_models/spec_severity_full.rds")

i_dis <- readRDS("../fit_models/indiv_disease_full.rds")
i_mort <- readRDS("../fit_models/indiv_mort_full.rds")
i_sev <- readRDS("../fit_models/indiv_severity_full.rds")


padding <- theme(plot.margin = unit(c(t=0.1, r=0.1, b=0, l=0), "cm"))

cp <-  egg::ggarrange(        
        (compare_posteriors(s_dis, i_dis) + padding + theme(axis.text.y=element_text(size=12)) + ggtitle("A) Disease presence")), 
        (compare_posteriors(s_mort, i_mort) + padding + theme(axis.text.y=element_blank()) + ggtitle("B) Mortality")), 
        (compare_posteriors(s_sev, i_sev) + padding + theme(axis.text.y=element_blank(), legend.position = "right") + ggtitle("C) Severity")),
         ncol=3, padding=0.5)

ggsave("../plots_tables/combined_species_individual_intervals.pdf", cp, width=11.5, height=3.5)

Human Case Fatality Rates

Host_order mean_cfr median_cfr min_cfr max_cfr
Rodentia 7.57 1.17 0 100
Chiroptera 62.90 80.92 0 100

CFR models

Comparing CFR for all viruses, regardless of host origin, and excluding bat viruses.

dat_cfr$Virus_name <- dat_cfr$Virus_ICTV
dat_cfr$Host_name <- dat_cfr$Host_Upham

# do bats and rodents differ in mean human CFR of inoculated viruses?

# model with no taxonomic / phylogenetic effects - rodents infected with lower human cfr viruses

if (!file.exists("../fit_models/cfr_base1.rds")) {

    cfr_base1 <- brm(CFR_avg/100 ~ Host_order+ 
                          (1|PaperID), 
        data=dat_cfr, family=zero_one_inflated_beta(),  
        prior=custom_prior,
        iter=4000, thin=1,
        control=list(adapt_delta=0.85, max_treedepth=12), cores=4, 
        save_pars=save_pars(all=TRUE))

    saveRDS(cfr_base1, "../fit_models/cfr_base1.rds")

} else { cfr_base1 <- readRDS("../fit_models/cfr_base1.rds")}

# cfr_base1
# pp_check(cfr_base1)
# cond_eff_plot(cfr_base1)
# param_tab(cfr_base1)

posterior <- as.matrix(cfr_base1)

color_scheme_set("teal")

ar1 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean", 
        point_size = 3.5, outer_size=1, inner_size=3) 
# ar1


# model with only viral taxonomic hierarcical effects

if (!file.exists("../fit_models/cfr_base2.rds")) {
   
    cfr_base2 <- brm(CFR_avg/100 ~ Host_order + 
                          (1|Virus_name) +
                          (1|Virus_ICTV) + 
                          (1|PaperID), 
        data=dat_cfr, family=zero_one_inflated_beta(),  
        prior=custom_prior,
        iter=5000, thin=1, cov_ranef = list(Virus_ICTV=virus_cov),
        control=list(adapt_delta=0.90, max_treedepth=12), cores=4, 
        save_pars=save_pars(all=TRUE))

    saveRDS(cfr_base2, "../fit_models/cfr_base2.rds")

} else { cfr_base2 <- readRDS("../fit_models/cfr_base2.rds")}

# cfr_base2

# pp_check(cfr_base2, ndraws = 200)
# cond_eff_plot(cfr_base2)
# param_tab(cfr_base2)

color_scheme_set("teal")
posterior <- as.matrix(cfr_base2)
ar2 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean", 
        point_size = 3.5, outer_size=1, inner_size=3)
# ar2


# model with no viral effects, but removing bat-reservoired viruses

bat_res_viruses <- vtraits$Virus_ICTV[vtraits$Reservoir%in%"Chiroptera"] %>% unique()
dat_cfr_nobatv <- dat_cfr[!dat_cfr$Virus_ICTV%in%bat_res_viruses,]

if (!file.exists("../fit_models/cfr_base3.rds")) {

    cfr_nobatv <- brm(CFR_avg/100 ~ Host_order + 
                          (1|PaperID), 
        data=dat_cfr_nobatv,  family=zero_one_inflated_beta(),  
        prior=custom_prior,
        iter=4000, thin=1,
        control=list(adapt_delta=0.85, max_treedepth=12), cores=4, 
        save_pars=save_pars(all=TRUE))

    saveRDS(cfr_nobatv, "../fit_models/cfr_base3.rds")

} else { cfr_nobatv <- readRDS("../fit_models/cfr_base3.rds")}

# cfr_nobatv

# pp_check(cfr_nobatv, ndraws=500)
# cond_eff_plot(cfr_nobatv)
# param_tab(cfr_nobatv)

posterior <- as.matrix(cfr_nobatv)
ar3 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean", 
        point_size = 3.5, outer_size=1, inner_size=3)
# ar3 


# # combo plot
# p1 <- plot_grid(  ar1 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nAll viruses"), 
#                   ar3 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nExcluding bat-reservoired viruses"), 
#                   ar2 + xlim(-3, 2) + scale_y_discrete(labels = "With virus effects\nAll viruses"),
#                   ncol=1, align = "v")
# p2 <- add_sub(p1, "Host order (Chiroptera) effect", x = 0, hjust = -3.5, size=10)
# ggdraw(p2)

# ggsave("../plots_tables/CFR_models_intervals.pdf", width=6, height=4)




### Models with only heterologous hosts

dat_cfr_het <- left_join(dat_cfr, dat %>% select(Host_Upham, Virus_ICTV, Reservoir_match, PaperID) %>% unique())
dat_cfr_het <- dat_cfr_het[dat_cfr_het$Reservoir_match==0,]


# model with no virus hierarchical effects

if (!file.exists("../fit_models/cfr_het1.rds")) {

    cfr_het1 <- brm(CFR_avg/100 ~ Host_order + (1|PaperID), 
        data=dat_cfr_het,  family=zero_one_inflated_beta(), 
        prior=custom_prior,
        iter=4000, thin=1,
        control=list(adapt_delta=0.85, max_treedepth=12), cores=4, 
        save_pars=save_pars(all=TRUE))

    saveRDS(cfr_het1, "../fit_models/cfr_het1.rds")

} else { cfr_het1 <- readRDS("../fit_models/cfr_het1.rds")}

# cfr_het1

# pp_check(cfr_het1, ndraws=300)
# cond_eff_plot(cfr_het1)
# param_tab(cfr_het1)

posterior <- as.matrix(cfr_het1)
ar4 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean", 
        point_size = 3.5, outer_size=1, inner_size=3)
# ar4 


# model with only viral taxonomic hierarcical effects

if (!file.exists("../fit_models/cfr_het2.rds")) {
   
    cfr_het2 <- brm(CFR_avg/100 ~ Host_order + 
                          (1|Virus_name) + 
                          (1|Virus_ICTV) + 
                          (1|PaperID),
        data=dat_cfr_het, family=zero_one_inflated_beta(),  
        prior=custom_prior,
        iter=5000, thin=1, cov_ranef = list(Virus_ICTV=virus_cov),
        control=list(adapt_delta=0.90, max_treedepth=16), cores=4, 
        save_pars=save_pars(all=TRUE))

    saveRDS(cfr_het2, "../fit_models/cfr_het2.rds")

# 0 div trans

} else { cfr_het2 <- readRDS("../fit_models/cfr_het2.rds")}

# pairs(cfr_het2)

# cond_eff_plot(cfr_het2)
# param_tab(cfr_het2)

color_scheme_set("teal")
posterior <- as.matrix(cfr_het2)
ar5 <- mcmc_intervals(posterior, regex_pars = "b_[^I]", prob=0.8, prob_outer=0.95, point_est="mean", 
        point_size = 3.5, outer_size=1, inner_size=3)
# ar5

# combo plot
p1 <- plot_grid(  ar1 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nAll viruses"), 
                  ar3 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nExcluding bat-reservoired viruses"), 
                  ar2 + xlim(-3, 2) + scale_y_discrete(labels = "With virus effects\nAll viruses"),
                  ar4 + xlim(-3, 2) + scale_y_discrete(labels = "No virus effects\nHeterologous hosts"),
                  ar5 + xlim(-3, 2) + scale_y_discrete(labels = "With virus effects\nHeterologous hosts"),
                  ncol=1, align = "v")
p2 <- add_sub(p1, "Host order (Chiroptera) effect", x = 0, hjust = -3.5, size=10)
ggdraw(p2)

ggsave("../plots_tables/CFR_models_intervals.pdf", width=6, height=6)

# dat_cfr %>% ungroup() %>% select(Virus_ICTV, CFR_avg) %>% unique() %>% arrange(CFR_avg) %>% View()